Computation of a STAI dataset with Field II using parameters of an L11-4v 128 element Verasonics Transducer and beamforming with USTB

This example shows how to load the data from a Field II simulation into USTB objects, and then beamform it with the USTB routines. This example uses the L11-4v 128 element Verasonics Transducer The Field II simulation program (field-ii.dk) should be in MATLAB's path.

This tutorial assumes familiarity with the contents of the 'CPWC simulation with the USTB built-in Fresnel simulator' tutorial. Please feel free to refer back to that for more details.

authors: Alfonso Rodriguez-Molares alfonso.r.molares@ntnu.no Ole Marius Hoel Rindal olemarius@olemarius.net

Last updated: 07.08.2017

Contents

Clear old workspace and close old plots

close all;

Basic Constants

Our first step is to define some basic constants for our imaging scenario - below, we set the speed of sound in the tissue, sampling frequency and sampling step size in time.

c0=1540;     % Speed of sound [m/s]
fs=100e6;    % Sampling frequency [Hz]
dt=1/fs;     % Sampling step [s]

field II initialisation

Next, we initialize the field II toolbox. Again, this only works if the Field II simulation program (field-ii.dk) is in MATLAB's path. We also pass our set constants to it.

field_init(0);
set_field('c',c0);              % Speed of sound [m/s]
set_field('fs',fs);             % Sampling frequency [Hz]
set_field('use_rectangles',1);  % use rectangular elements
      *------------------------------------------------------------*
      *                                                            *
      *                      F I E L D   I I                       *
      *                                                            *
      *              Simulator for ultrasound systems              *
      *                                                            *
      *             Copyright by Joergen Arendt Jensen             *
      *    Version 3.30, April 5, 2021 (Matlab 2021a version)      *
      *                  Web-site: field-ii.dk                     *
      *                                                            *
      *     This is citationware. Note the terms and conditions    *
      *     for use on the web-site at:                            *
      *               field-ii.dk/?copyright.html                  *
      *  It is illegal to use this program, if the rules in the    *
      *  copyright statement is not followed.                      *
      *------------------------------------------------------------*
Warning:  Remember to set all pulses in apertures for the new sampling frequency

Transducer definition L11-4v, 128-element linear array transducer

Our next step is to define the ultrasound transducer array we are using. For this experiment, we shall use the L11-4v 128 element Verasonics Transducer and set our parameters to match it.

probe = uff.linear_array();
f0                      = 5.1333e+06;      % Transducer center frequency [Hz]
lambda                  = c0/f0;           % Wavelength [m]
probe.element_height    = 5e-3;            % Height of element [m]
probe.pitch             = 0.300e-3;        % probe.pitch [m]
kerf                    = 0.03e-03;        % gap between elements [m]
probe.element_width     = probe.pitch-kerf;% Width of element [m]
lens_el                 = 20e-3;           % position of the elevation focus
probe.N                 = 128;             % Number of elements
pulse_duration          = 2.5;             % pulse duration [cycles]

Pulse definition

We then define the pulse-echo signal which is done here using the fresnel simulator's pulse structure. We could also use 'Field II' for a more accurate model.

pulse = uff.pulse();
puse.center_frequency = f0;
pulse.fractional_bandwidth = 0.65;        % probe bandwidth [1]
t0 = (-1/pulse.fractional_bandwidth/f0): dt : (1/pulse.fractional_bandwidth/f0);
impulse_response = gauspuls(t0, f0, pulse.fractional_bandwidth);
impulse_response = impulse_response-mean(impulse_response); % To get rid of DC

te = (-pulse_duration/2/f0): dt : (pulse_duration/2/f0);
excitation = square(2*pi*f0*te+pi/2);
one_way_ir = conv(impulse_response,excitation);
two_way_ir = conv(one_way_ir,impulse_response);
lag = length(two_way_ir)/2+1;

% We display the pulse to check that the lag estimation is on place
% (and that the pulse is symmetric)

figure;
plot((0:(length(two_way_ir)-1))*dt -lag*dt,two_way_ir); hold on; grid on; axis tight
plot((0:(length(two_way_ir)-1))*dt -lag*dt,abs(hilbert(two_way_ir)),'r')
plot([0 0],[min(two_way_ir) max(two_way_ir)],'g');
legend('2-ways pulse','Envelope','Estimated lag');
title('2-ways impulse response Field II');

Aperture Objects

Next, we define the the mesh geometry with the help of Field II's xdc_linear_array function.

noSubAz=round(probe.element_width/(lambda/8));        % number of subelements in the azimuth direction
noSubEl=round(probe.element_height/(lambda/8));       % number of subelements in the elevation direction
Th = xdc_linear_array (probe.N, probe.element_width, probe.element_height, kerf, noSubAz, noSubEl, [0 0 Inf]);
Rh = xdc_linear_array (probe.N, probe.element_width, probe.element_height, kerf, noSubAz, noSubEl, [0 0 Inf]);

% We also set the excitation, impulse response and baffle as below:
xdc_excitation (Th, excitation);
xdc_impulse (Th, impulse_response);
xdc_baffle(Th, 0);
xdc_center_focus(Th,[0 0 0]);
xdc_impulse (Rh, impulse_response);
xdc_baffle(Rh, 0);
xdc_center_focus(Rh,[0 0 0]);

Phantom

In our next step, we define our phantom. Here, our phantom is a single point scatterer.

g2=[[0.1e-3 0.2e-3 0.5e-3 1e-3 2e-3].' zeros(5,1) linspace(22e-3,18e-3,5).'];
g1=[zeros(5,1)              zeros(5,1) linspace(18e-3,22e-3,5).'];             % list with the scatterers coordinates [m]
sca=[g1; g2];
amp=ones(10,1);                                                           % list with the scatterers amplitudes

Output data

We define the variables to store our output data

cropat=round(1.1*2*sqrt((max(sca(:,1))-min(probe.x))^2+max(sca(:,3))^2)/c0/dt);   % maximum time sample, samples after this will be dumped
STA=zeros(cropat,probe.N,probe.N);    % impulse response channel data

Compute STA signals

Now, we finally reach the stage where we generate a STA (Synthetic Transmit Aperture) dataset with the help of Field II.

disp('Field II: Computing STA dataset');
wb = waitbar(0, 'Field II: Computing STA dataset');
for n=1:probe.N
    waitbar(n/probe.N, wb);

    % transmit aperture
    xdc_apodization(Th, 0, [zeros(1,n-1) 1 zeros(1,probe.N-n)]);
    xdc_focus_times(Th, 0, zeros(1,probe.N));

    % receive aperture
    xdc_apodization(Rh, 0, ones(1,probe.N));
    xdc_focus_times(Rh, 0, zeros(1,probe.N));

    % do calculation
    [v,t]=calc_scat_multi(Th, Rh, sca, amp);

    % build the dataset
    STA(1:size(v,1),:,n)=v;

    % Sequence generation
    seq(n)=uff.wave();
    seq(n).probe=probe;
    seq(n).source.xyz=[probe.x(n) probe.y(n) probe.z(n)];
    seq(n).sound_speed=c0;
    seq(n).delay = probe.r(n)/c0-lag*dt+t; % t0 and center of pulse compensation
end
close(wb);
Index exceeds the number of array elements (1).

Error in STAI_L11_resolution_phantom (line 127)
disp('Field II: Computing STA dataset');

Channel Data

In this part of the code, we creat a uff data structure to specifically store the captured ultrasound channel data.

channel_data = uff.channel_data();
channel_data.sampling_frequency = fs;
channel_data.sound_speed = c0;
channel_data.initial_time = 0;
channel_data.pulse = pulse;
channel_data.probe = probe;
channel_data.sequence = seq;
channel_data.data = STA./max(STA(:));

Scan

The scan area is defines as a collection of pixels spanning our region of interest. For our example here, we use the linear_scan structure, which is defined with two components: the lateral range and the depth range. scan too has a useful plot method it can call.

SCAN

sca=uff.linear_scan('x_axis',linspace(-2e-3,4e-3,256).', 'z_axis', linspace(17e-3,23e-3,256).');

Pipeline

With channel_data and a scan we have all we need to produce an ultrasound image. We now use a USTB structure pipeline, that takes an apodization structure in addition to the channel_data and scan.

pipe=pipeline();
pipe.channel_data=channel_data;
pipe.scan=sca;

pipe.receive_apodization.window=uff.window.boxcar;
pipe.receive_apodization.f_number=1.7;
pipe.transmit_apodization.window=uff.window.boxcar;
pipe.transmit_apodization.f_number=1.7;

The pipeline structure allows you to implement different beamformers by combination of multiple built-in processes. By changing the process chain other beamforming sequences can be implemented. It returns yet another UFF structure: beamformed_data.

To achieve the goal of this example, we use delay-and-sum (implemented in the das_mex() process) as well as coherent compounding.

b_data=pipe.go({midprocess.das() postprocess.coherent_compounding()});

% Display images
b_data.plot()

Save UFF dataset

Finally, we save the data into a UFF file.

%channel_data.write([ustb_path(),'/data/FieldII_PSF_simulation_v2.uff'],'channel_data');